

##############################################　Survey Experiment on Diversionary War code by Shingo Hamanaka 20221104 ###############################################
IsraelList <- read.csv("~/IsraelList.csv")
### reading datafile

library(ggplot2)
library(list)
library(tidyverse)
#################### missing variables ##########################################
data1 <-IsraelList[!(is.na(IsraelList$listoutcome) |  is.na(IsraelList$Thermo_Bibi) | is.na(IsraelList$Patriotism) | is.na(IsraelList$experiment) | is.na(IsraelList$support_Bibi) | is.na(IsraelList$Gender) | is.na(IsraelList$Age_Group) | is.na(IsraelList$Education) | is.na(IsraelList$Income) | is.na(IsraelList$Ideology)),]


################### full model #############################

diversion <- ictreg.joint(formula = listoutcome ~ Thermo_Bibi + Patriotism + Ideology + Gender + Income + Education + Age_Group, data=data1, treat="experiment", outcome="support_Bibi", J=3, constrained=TRUE, outcome.reg="linear", maxIter=1000)

summary(diversion)

Bibipred <- predict.ictreg.joint(diversion, se.fit = TRUE, interval = "confidence",level = 0.95, avg = TRUE,
                                     sensitive.value = "both",
                                     sensitive.diff = TRUE, return.draws = TRUE,
                                     predict.sensitive = TRUE)
Bibipred$fit

Bibipred$fitsens

as_tibble(Bibipred$fit) %>% 
  mutate(term = c("Diversion disagree", "Diversion agree", "Difference") %>% 
           fct_inorder()) %>% 
  ggplot() + 
  geom_pointrange(aes(x = term, y = fit, ymin = lwr, ymax = upr)) + labs(title="Full Model") +
  geom_hline(aes(yintercept = 0), linetype = "dashed")

################### simple model #############################
simplemodel <- ictreg.joint(formula = listoutcome ~ Gender + Income + Education + Age_Group, data=data1, treat="experiment", outcome="support_Bibi", J=3, constrained=TRUE, outcome.reg="linear", maxIter=1000)

summary(simplemodel)

simpleBibipred <- predict.ictreg.joint(simplemodel, se.fit = TRUE, interval = "confidence",level = 0.95, avg = TRUE,
                                       sensitive.value = "both",
                                       sensitive.diff = TRUE, return.draws = TRUE,
                                       predict.sensitive = TRUE)
simpleBibipred$fit

simpleBibipred$fitsens

as_tibble(simpleBibipred$fit) %>% 
  mutate(term = c("Diversion disagree", "Diversion agree", "Difference") %>% 
           fct_inorder()) %>% 
  ggplot() + 
  geom_pointrange(aes(x = term, y = fit, ymin = lwr, ymax = upr)) + labs(title="Simple Model") +
  geom_hline(aes(yintercept = 0), linetype = "dashed")
